package gene_loci;


require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(copy_ordered_gene_loci_by_chromosome);

use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use strict;
use Carp;
use warnings;
use chromosomes;
use constant GENEID =>0;
use constant CHRMSM =>1;
use constant START_ =>2;
use constant END___ =>3;
use constant STRAND =>4;
use constant CENTRE =>5;
use constant ORDER_ =>5;

#_________________________________________________________________________________________

#	merge_overlapping_transcripts

#		merge the gene loci of different transcripts from the same gene together


#_________________________________________________________________________________________
sub merge_overlapping_transcripts(\@\@)
{
	# array of array of $gene_id, $chromosome, $start, $end, $strand
	my ($gene_loci, $transcript_strand_errors) = @_;

	croak "Error\n\tNo gene at loci\n" unless @$gene_loci;

	#
	# initialize with first loci
	#
	my ($gene_id, $chromosome, $start, $end, $strand) =
									@{$gene_loci->[0]}[GENEID, CHRMSM, START_, END___, STRAND];

	#
	# go through all other loci
	#
	for my $i(1 .. (@$gene_loci - 1))
	{
		#
		# complain if different chromosomes
		#
		if ($chromosome ne $gene_loci->[$i][CHRMSM])
		{

			#
			# Handle 'Un' and 'random' chromosome. These are unplaced reads.
			#	The loci are completely useless and should be superceded
			#	by anything else
			#
			if ($chromosome eq 'Un')
			{
				# supercede old unknown
				($chromosome, $start, $end, $strand) =
										@{$gene_loci->[$i]}[CHRMSM, START_, END___, STRAND];
				next;
			}
			elsif ($chromosome =~/(.*)_random/)
			{
				#
				# if currently saved is random and encounter non-random
				# have non-random over-write random
				my $base_chromosome = $1;
				if ($gene_loci->[$i][CHRMSM] eq $base_chromosome)
				{
					($chromosome, $start, $end, $strand) =
											@{$gene_loci->[$i]}[CHRMSM, START_, END___, STRAND];
					next;
				}
			}
			elsif ($gene_loci->[$i][CHRMSM] eq 'Un')
			{
				# ignore new unknown
				next;
			}
			elsif ($gene_loci->[$i][CHRMSM] =~/(.*)_random/)
			{
				#
				# if currently saved is non-random and encounter random
				# ignore that
				my $base_chromosome = $1;
				next if ($chromosome eq $base_chromosome);
			}

			warn (('!' x80).
				 "\n\nWARNING:!!!!\n\n".
				 "\tTranscripts for [$gene_id] have different chromosomes ".
				"($chromosome vs. $gene_loci->[$i][CHRMSM])\n".
				 ('!' x80)."\n");
			next;
		}

		push(@$transcript_strand_errors,  $gene_id)
			 unless $strand	== $gene_loci->[$i][STRAND];
		$start = $gene_loci->[$i][START_] if $start > $gene_loci->[$i][START_];
		$end = $gene_loci->[$i][END___] if $end < $gene_loci->[$i][END___];
	}

	return [$gene_id, $chromosome, $start, $end, $strand, ($start + $end) / 2];
}

#_________________________________________________________________________________________

#	order_gene_loci_by_chromosome

#			merge transcripts and sort into chromosome buckets

#_________________________________________________________________________________________
sub copy_ordered_gene_loci_by_chromosome($\@\@\@$)
{
	my ($dbh, $gene_loci, $transcript_strand_errors, $transcript_length_errors, $db_origin) = @_;

	print STDERR "\t\tSorting ".scalar(@$gene_loci)." transcripts...\n";
	#
	#	group transcripts together by gene_id
	#
	my %transcripts_by_gene_id;
	for my $transcript_loci(@$gene_loci)
	{
		push(@{$transcripts_by_gene_id{$transcript_loci->[GENEID]}}, $transcript_loci);
	}

	#
	#	merge transcripts with different loci for the same gene_id
	#
	my %chromosomes;
	my %gene_loci_by_chromosome;
	for my $loci( keys %transcripts_by_gene_id)
	{

		my $locus = merge_overlapping_transcripts(@{$transcripts_by_gene_id{$loci}},
													@$transcript_strand_errors);
		# save by chromosome
		push(@{$gene_loci_by_chromosome{$locus->[CHRMSM]}}, $locus);
		#save list of all chromosomes
		++$chromosomes{$locus->[CHRMSM]};
	}

	#
	#	sort chromosomes
	#
	my @sorted_chromosomes = sort_chromosomes(keys %chromosomes);


	#
	#	sort within chromosomes
	#
	for my $chromosome(@sorted_chromosomes)
	{
		#
		#	sort by centre point of gene loci
		#	then by start
		#	then by end
		#	then by geneid if all else fails!
		#
		@{$gene_loci_by_chromosome{$chromosome}} = sort{
														$a->[CENTRE] <=> $b->[CENTRE] or
														$a->[START_] <=> $b->[START_] or
														$a->[END___] <=> $b->[END___] or
														$a->[GENEID] cmp $b->[GENEID]
													}
												@{$gene_loci_by_chromosome{$chromosome}};
	}

	my @all_sorted_gene_loci;
	my $gene_order = 0;
	for my $chromosome(@sorted_chromosomes)
	{
		for my $gene_locus(@{$gene_loci_by_chromosome{$chromosome}})
		{
			# save list of genes which span too much of the chromosome
			if ($gene_locus->[END___] - $gene_locus->[START_] > 1000000)
			{
				# if too long: Save
				push(@$transcript_length_errors,
							[$gene_locus->[GENEID], 							# gene_id
							($gene_locus->[END___] - $gene_locus->[START_])]);	# length
			}

			my $data = "$gene_locus->[GENEID]\t".
						"$gene_locus->[CHRMSM]\t".
						"$gene_locus->[START_]\t".
						"$gene_locus->[END___]\t".
						($gene_locus->[STRAND] == 1 ? 1 : 0)."\t".
						(++$gene_order)."\t".
						"$db_origin\n";
			($$dbh)->pg_putline($data);
		}
	}
	print STDERR "\t\t$gene_order gene loci copied to panda...\n";
}






1;
